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Abstract. Solving analytic systems using inversion can be implemented 
in a variety of ways. One method is to use Lagrange inversion and varia- 
tions. Here we present a different approach, based on dual vector fields. 

For a function analytic in a neighborhood of the origin in the complex 
plane, we associate a vector field and its dual, an operator version of 
Fourier transform. The construction extends naturally to functions of 
several variables. 

We illustrate with various examples and present an efficient algorithm 
readily implemented as a symbolic procedure in Maple while suitable as 
well for numerical computations using languages such as C or Java. 



1 Introduction 

We introduce the operator calculus necessary to present our approach to (local) 
inversion of analytic functions. It is important to note that this is different from 
Lagrange inversion and is based on the flow of a vector field associated to a given 
function. It appears to be theoretically appealing as well as computationally 
effective. 

Acting on polynomials in x, define the operators 

D = — and X = multiplication by x. 
dx 

They satisfy commutation relations [D, X] = I, where /, the identity opera- 
tor, commutes with both D and X. Abstractly, the Heisenberg-Weyl algebra is 
the associative algebra generated by operators {A, B, C} satisfying [A, B] — C, 
[A, C] — [B,C] = 0. The standard HW algebra is the one generated by the 
realization A = D, B = X, C = I. An Appell system is a system of polynomials 
{y n (x)} n >o that is a basis for a representation of the standard HW algebra with 
the following properties: 

1. y n is of degree n in x; 

2. Dy n = ny n -i- 

In several variables, x = (x±, . . . , xn), with multi-indices n = (m, . . . , njsr), the 
corresponding monomials are 
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Denote the partial derivative operators by D{ = — — and the corresponding 

OXi 

multiplication operators by Xi. Then [Dj,Xi] = SijI. An Appell system is a 
system of polynomials {y n } in the variables x such that 

1. the top degree term of y n is a constant multiple of x n ; 

2- Diy n — 7^i^/n-ei, where has all components zero except for 1 in the i th 
position. 

G.-C. Rota [H| is well-known for his umbral calculus development of special 
polynomial sequences, called basic sequences. From our perspective, these are 
"canonical polynomial systems" in the sense that they provide polynomial rep- 
resentations of the Heisenberg-Weyl algebra, in realizations different from the 
standard one. Our idea |2llj is to illustrate explicitly the role of vector fields and 
their duals, using operator calculus methods for working with the latter (in our 
volumes — this viewpoint is prefigured in |3J). 

The main feature of our approach is that the action of the vector field may 
be readily calculated while the action of the dual vector field on exponentials is 
identical to that of the vector field. Then we note that acting iteratively with 
a vector field on polynomials involves the complexity of the coefficients, while 
acting iteratively with the dual vector field always produces polynomials from 
polynomials. So we can switch to the dual vector field for calculations. 

Specifically, fix a neighborhood of in C. Take an analytic function V(z) 
defined there, normalized to V(0) = 0, V'(0) = 1. Denote W(z) = l/V'(z) and 
U(v) the inverse function, i.e., V(U(v)) = v, U(V(z)) = z. Then V(D) is defined 
by power series as an operator on polynomials in x and [V(Z)),A] = V'(D) 
so that [V(D),XW(D)] = I. In other words, V = V(D) and Y = XW(D) 
generate a representation of the HW algebra on polynomials in x. The basis for 
the representation is y n {x) = Y n l, i.e., Y is a raising operator. And Vy n — 
ny n -i so that V is the corresponding lowering operator. The {y n }n>o form a 
system of canonical polynomials or generalized Appell system. The operator of 
multiplication by x is given by X = YV'(D) — YU' (V)^ 1 , which is a recursion 
operator for the system. 

We identify vector fields with first-order partial differential operators. Con- 
sider a variable A with corresponding partial differential operator 8a- Given V 
as above, let Y be the vector field Y = W(A)dA- Then we observe the following 
identities 

Y e Ax = xW(A) e Ax = xW{D) e Ax 

as any operator function of D acts as a multiplication operator on e Ax . The 
important property of these equalities is that Y and Y commute, as they involve 
independent variables. So we may iterate to get 



eMtY)e Ax = eMtY)e Ax . 



(1) 
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On the other hand, we can solve for the left-hand side of this equation using 
the method of characteristics. Namely, if we solve 

A = W(A) (2) 

with initial condition A(0) = A, then for any smooth function /, 

e*f{A) = f(A(t)). 

Thus 

exp(tY)e Ax =e xA{t) . 
To solve equation ©, multiply both sides by V'(A) and observe that we get 

V'(A)A = j t V{A{t)) = l. 

Integrating yields 

V(A(t)) = t + V(A) or A(t) = U(t + V(A)). 
Or, writing v for t, we have 

exp(vY)e Ax = e wU ^+ v ^K (3) 

We can set A — to get 

exp{vY)l =e l(/w 

on the one hand while 



e vY, 



1 = TVn{x). 
^-^ TV. 

71=0 

In summary, we have the expansion of the exponential of the inverse function 



n=0 



or 



E^r(^)) m = E^^)- ( 4 ) 



mi * — ' n 

m—0 n—Q 



This yields an alternative approach to inversion of the function V(z) rather 
than using Lagrange's formula. We see that the coefficient of x m /ml yields the 
expansion of (U(v)) m . In particular, U(v) itself is given by the coefficient of x 
on the right-hand side. 
Specifically, we have: 

Theorem 1. The coefficient of x m /m\ in Y n l is equal to Y n A m \ A _ Q , each 
giving the coefficient ofv n /n\ in the expansion ofU(v) m . 
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Proof. Expand both sides of equation Q, using v for t, in powers of x and v, 
and let A = 0: 

n! ^— ' ml ,A =° ^— ' ra! 

n— m— n=0 

and compare with equation Q. 

The same idea works in several variables. 

We have V(z) = (Vi(zi, . . . , zjy), •••■> Vtf(zi, . . . , z/v)) analytic in a neighbor- 
hood of in C N . Denote the Jacobian matrix [ j by V' and its inverse by 

\dzjj 

W. The variables 

JV 

Yi = J2 x kW ki (D) 
fc=i 

commute and act as raising operators for generating the basis y n (x). Namely 
YiVn = 2/n+e;- And Vi = Vi(D), D = (Di, . . . , Dn) are lowering operators: 
Viy n = rii y n - ei . 

Denote ^2 i aibi by a ■ b. With variables Ai and corresponding partials di, 
define the vector fields 

fc 

For a vector field Y — ^\ we have the identities 

Y e A - x = x ■ W{A) e A ' x = x ■ W(D) e A ' x . 
The method of characteristics applies as in one variable and as in equation J3J) 

eMvY)e Ax = e *-u(v+v { A))_ 

Thus, we have the expansion 

exp(a; • U («)) = X, Z/n(x). (5) 

n 

In particular, the fc th component, of the inverse function is given by the 
coefficient of Xk in the above expansion. 

An important feature of our approach is that to get an expansion to a given 
order requires knowledge of the expansion of W just to that order. The reason 
is that when iterating xW(D), at step n it is acting on a polynomial of degree 
n — 1, so all terms of the expansion of W(D) of order n or higher would yield 
zero acting on y n -i- This allows for streamlined computations. 

For polynomial systems V, V' will have polynomial entries, and W will be 
rational in z. Hence raising operators will be rational functions of D, linear in x. 
Thus the coefficients of the expansion of the entries of W would be computed 
by finite-step recurrences. 

Remark 1. Note that to solve V{z) = v for z near zq, with V(zq) — vq, apply 
the method to V\{z) = V(z + z ) — v , so that Vi(0) = 0. The inverse is U\{v) = 
U(v + v ) — Zq. Then U(v) = z + U\{v — Vq). 
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2 One-variable Case 

In this section we focus on the one-variable case. We illustrate the method with 
examples, and then present an algorithm suitable for symbolic computation. 

Example 1. In one variable, solving a cubic is interesting as the expansion of W 
can be expressed in terms of Chebyshev polynomials. 
Let V = z 3 /3 - az 2 + z. Then V 



1 — 2az + z 

7L — U 

where U n are Chebyshev polynomials of the second kind. 

Specializing a provides interesting cases. For example, let a — cos(7r/4), or 
V = z 3 /3 — z 2 /y/2 + z. Then the coefficients in the expansion of W are periodic 
with period 8 and, in fact, 

W 1 + ^ + 



1 + z 4 

The coefficient of x in the polynomials y n yield the coefficients in the expansion 
of the inverse U . Here are some polynomials starting with j/o = 1, j/i = x: 

y 2 =x 2 +xV2, y 3 = x 3 + 3x 2 V2 + Ax, 

?/4 = x 4 + 6 x 3 V2 + 22 x 2 + 10 x V2, 

y 5 = x 5 + 10 x 4 V2 + 70 x 3 + 90 x 2 V2 + 40 x, 

y 6 = x 6 + 15 x 5 V2 + 170 x* + 420 x 3 V2 + 700 x 2 - 140 x V2. 

This gives to order 6: 



This expansion will give approximate solutions to 

z 3 /3~z 2 /V2 + z-v = 

for v near 0. 

Example 2. Inversion of the Chebyshev polynomial Tz{z) = Az 3 — 3z can be used 
as the basis for solving general cubic equations (0). 

To get started we have, with V(z) — Az 3 — 3z, 



W(z) = — L_j = — V A n z 2n . 

y ' 3 1-Az 2 3 ^ 

So Vl = (-l/3)a;, 2/2 - (l/9)a; 2 , y 3 = {-l/27)(x 3 + 8a:), etc. We find 
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In this case, we can find the expansion analytically. To solve T 3 (z) — v, write 

T 3 (cos6>) = cos(3<9) = v. 

Invert to get, for integer k, 9 = (l/3)(27rfc ± arccosw), with arccos denoting the 
principal branch. Then 

z = cos((l/3)(27rA: ± arccos i;)). 

We want a branch with v = corresponding to z = 0. With arccos = n/2, we 
want the argument of the cosine to be n/2 + ttI, for some integer I. This yields 
1 21 + 1 

the condition - = — -. Taking I = 0, we get k = 1, with the minus sign. 

3 4fc ± 1 

Namely, 

U(v) = cos((l/3)(27r — arccos v)). 

Using hypergeometric functions (see next example) and rewriting, we find the 
form 

'3n\ /4\" v 2n+1 



u( V ) = 4 e ( 



. n J V 27 y 2n+l 

n=0 x ' x ' 

If we generate the polynomials y n , we can find the expansion of U(v) m to 
any order. 

Example 3. A similar approach is interesting for the Chebyshev polynomial 
T n (z). 

F(v) — cos( A (//± arccos v)) satisfies the hypergeometric differential equation 

(1 - v 2 ) F" - v F' + X 2 F = 

which can be written in the form 

[{vD v f - D 2 V ]F = A 2 F 

with here D v denoting d/dv. For integer A, this is the differential equation for the 
corresponding Chebyshev polynomial. In general, these are Chebyshev functions. 
As noted above, for F(0) = 0, we take fi — 2irk, and, as above, we require 

21 + 1 
4k ± 1 ' 

With F'(0) = ±A, we have the solution 



F(v) = ±\v 2 F 1 



/ 1 + A 1- A 

2 ' 2 
3 

V 2 
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2.1 Using Maple 

For symbolic computation using Maple, one can use the Ore_Algebra package. 

1. First fix the degree of approximation. Expand W as a polynomial to that 
degree. 

2. Declare the Ore algebra with one variable, x, and one derivative, D. 

3. Define the operator xW(D) in the algebra. 

4. Iterate starting with yo = 1 using the applyopr command. 

5. Extract the coefficient of x m jm\ to get the expansion of U(v) m . 



3 Algorithm as a Matrix Computation 



Here is a matrix approach that can be implemented numerically 
Fix the order of approximation n. Cut off the expansion 



W(z) = WQ + W\Z + W2W 2 + 



w k z 



at w n z n . 

Let the matrix 



w 2 Wi 



W = 





w 



Wn-l W n - 2 W„_ 3 
\ W n W„_i W n -2 



0\ 




• w 

■ WiJ 



Define the auxiliary diagonal matrices 



Q 



Note that QP = M. 



(\\ ... 0\ 
2! ... 



\0 



M 
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(i/r(\) o ... 
o i/r(2)... 



/io 

2 



V o o 







\0 ... n) 

\ 





Denoting yk(x) = c ^ x3 , we have the recursion 



r (fe+i) (fe+i) 

Ci , to 



, . . . , v, n 



,...,c«]P^Q. 



The condition U(0) = gives j/o = 1- Then y\ = XW(D)y yields y\ 

(k) 

We see that c = for k > 0. We iterate as follows: 



WqX. 
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1. Start with w times the unit vector [1,0,..., 0] of length n. 

2. Multiply by W. 

3. Iterate, multiplying on the right by MW at each step. 

4. Finally, multiply on the right by Q. 

The top row will give the coefficients of the expansion of U(v) to order n. 

4 Higher-order Example 

Here is a simple 2x2 system for illustration. 

Vi = z\ + z%/2, V 2 — z 2 - Z\z 2 . 

So 

V'=( l / 2 ) and W=- l —( l - Z ^-^\. 

\-z 2 1 - ZiJ 1 - zi + z\ \ z 2 1 J 

The raising operators are 

Yi = (a;i(l-Di))+a:2l>2) (1 — Di + D 2 )~ l , 
Y 2 = {- Xl D 2 + x 2 ) (1 - £>! + D 2 2 )-\ 

OO 

Expanding (1 - D t + Df}- 1 = £ (D 1 - D\) n yields, with y 00 = 1, 

Uoi=x 2 , yw = x 1 , 
yo2 = x\-xi, yn=x 2 +x 1 x 2 , y 2 Q = x\. 

Thus 

exp(x • U(v)) = 1 + x\Vi + x 2 v 2 

+ (x 2 + xix 2 )viv 2 + {x\ -n)y + a;?yH , 

so 

Ui(v) =vi - vf/2 H , U 2 (v) =v 2 + viv 2 H . 

5 Another Matrix Approach 

For any given order n, the polynomials of degree n are an invariant subspace for 
the operator Y up until the last step. We can formulate an alternative matrix 
computation as follows. Let D and X denote the matrices of the operators of 
differentiation and multiplication by x respectively on polynomials of degree less 
than or equal to n. The space is invariant under differentiation, and we cut off 
multiplication by x to be zero on x n . We get 

Dij = iSi + ij and X t j = 
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with the first row of X all zeros. We then compute the matrix X times W(D), 
where W(D) is computed as a matrix polynomial by substituting in W{z) up to 
order n. Then Y has a matrix representation, Y = XW(D), on the space and we 
iterate multiplying by Y acting on the unit vector ei . These give the coefficients 
of the polynomials y n . 

In several variables, one constructs matrices for Dj and Xi using Kronecker 
products of D and X with the identity. For example, 

Dj = J®I®«--®.D®/'--®J 

with D in the j th spot. Similarly for Xi. Then one has explicit matrix represen- 
tations for the dual vector fields and the polynomials can be found accordingly. 

This approach is explicit, but seems to much slower than using the built-in 
Ore_algebra package. 

6 Worksheets 

> withCl inalg) : 

> withCOre_a1gebra) : 

One Variable 

_> 

> n:=10; 

n := 10 

_> unassignC'z' .'y'.'Y'.'x'.'V'.'W'): 

> V:=z-zA2/2; 

2 

> W:=diffCV,z)A(-l) ; 

> W:=convertCtay1or(W,z=0,n) .polynom) ; 

1 - z 

... 2 3 4 5 6 7 8 9 

W:= 1+z+z + z + z + z + z + z + z + z 

> ############ ORE ALGEBRA STARTS HERE ##################### 

> A:=diff_a1gebraC[z.x]) ; 

A := Ore_algebra 

> YY:=x*W; 

2345678 9. 
YY:= x (1 + z + z + z + z + z + z + z + z + Z ) 

> y[0]:=l: 

> for i from 1 to n-1 do y [i] :=s"impl i fy(applyoprCYY,y[i-l] , A)) od; 

y 1 := x 

2 

y 2 := x + x 

y 3 := x 3 + 3 x 2 + 3 x 

4 3 2 

y 4 :=x + 6x +15x + 15 x 

y 5 := x 5 + 10 x 4 + 45 x 3 + 105 x 2 + 105 x 
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Several Variables 

Example in two variables from an analytic function 

> VI : =eval c (Re(expand (subs (z=zc , f ) ) ) ) ; V2 : =eval c (Im(expand (subs (z=zc , f ) ) ) ) ; 

VI := zl-Azl 2 + 4zlz2+ 4 z2 2 

V2:= z2- 2zl 2 - 8 zl z2+ 2 z2 2 

> Jac : = j acobi an ( [VI , V2] , [zl , z2] ) : 

> W : =eval m( Jac a (- 1) ) :adj( Jac) .factor (det( Jac)) ; 

> WMat : =map(mtayl or , W, [zl=0 , z2=0] , n) ; 



1 - 8 zl + 4 z2 -4 zl - 8 z2 
4zl+8z2 1 - 8 zl + 4 z2 



1 - 16 zJ + 8z2+ 80 zl 2 + 80 z2 2 



WMat- [[1 + 8zl- 4z2 + 48 zl 2 - 48 z2 2 - 128 zl z2+ 128 zl 3 - 384 zl z2 2 

- I\\lz2zl 2 + 704 z2 3 , -4 zl- 8z2- 64 zl 2 - 96 zl z2+ 64 z2 2 - 704 zl 3 

- 384 z2zl 2 + 2112 zl z2 2 + 128 z2 3 ], [4 zl + 8 z2+ 64 zl 2 + % zl z2- 64 z2 2 
+ 704 zl 3 + 384 z2zl 2 - 2112 zl z2 2 - 128 z2 3 , 1 + 8 zl - 4z2+ 48 zl 2 - 48 z2 2 

- 128 zl z2+ 128 zl 3 - 384 zl z2 2 - l\\lz2zl 2 + 704 z2 3 ]] 

[> ############ ORE ALGEBRA STARTS HERE ##################### 

> A:=diff_algebra([zl,xl] , [z2,x2]) ; 

A:= Ore_algebra 

> for ix to N do YY[ix] :=simplify(add(x| |k*WMat[k,ix] ,k=l. .N)) od; 

YY 1 := xl+8xlzl-4xl z2+ 48x1 zl 2 -48xl z2 2 - 128 x1 zlz2+ 128x1 zl 3 

- 384x1 zl z2 2 - 2112x1 z2zl 2 + 704 xl z2 3 + 4x2 zl + 8x2 z2+ 64 x2zl 2 
+ 96 x2zl z2- 64 x2z2 2 + 704 x2zl 3 + 384 x2z2zl 2 - 2112x2 zl z2 2 

- 128 x2z2 3 

YY 2 := -4 xl zl - 8 xl z2- 64 xl zl 2 - 96 xl zl z2+ 64 xl z2 2 - 704 xl zl 3 

- 384x1 z2zl 2 + 2112x1 zl z2 2 + 128 x1 z2 3 + x2+ 8x2zl-4x2z2+ 48 x2zl 2 

- 48 x2z2 2 - 128 x2zl z2+ 128 x2zl 3 - 384 x2zl z2 2 - 2112 x2z2zl 2 
+ 704 x2 z2 3 

> y[0,0]:=l: 

> for i from to n-1 do j:=0; if(i>0) then 
y[i,0]:=simplify(applyopr(YY[l],y[i-l,0],A)) fi; for j from 1 to n-1 do 
y[i,j]:=simplify(applyopr(YY[2],y[i,j-l],A));print("y["||iM","||jM"]=", 
[i,j]) od; od; 

j:=0 
"y[0,l]=",x2 
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"y[0,l]=", x2 
"y[0,2]=", x2 2 - 4 x2- 8 xl 
"y[0,3]=", x2 3 - 12 x2 2 - 24 x2xl + 192 xl - 144 x2 
j:= 

"y[l,l]=", x2xl- 4x1+ 8x2 
"y[l,2]=", x2 2 xl - 12 x2xl + 16 x2 2 - 144 xl - 192 x2 - 8 xl 2 

"y[13]=", 10560x1 - 1920x2- 24 x2 2 xl + 24 x2 3 + x2 3 xl - 24 x2 xl 2 - 720x2x1 

- 672 x2 2 + 288 xl 2 

j:= 

"y[2,l]=", x2xl 2 + 24 x2 xl + 4 x2 2 - 8 xl 2 - 192 xl + 144 x2 
"y[2,2]=", -1920x1 - 10560x2+ 40 x2 2 xl - 8 xl 3 + 4 x2 3 + x2 2 xl 2 - 20 x2xl 2 

- 960 x2 xl + 400 x2 2 - 320 xl 2 

"y[2,3]=", 497664x1 + 145152x2+ x2 3 xl 2 - 24 x2 xl 3 - 2496 x2 2 xl + 384 xl 3 
+ 768 x2 3 - 36 x2 2 xl 2 + 56 x2 3 xl - 1392 x2 xl 2 + 4 x2 4 - 13440 x2 xl 

- 43200 x2 2 + 30720x1 2 

j:= 

"y[3,U=", x2xl 3 + 48 x2xl 2 + 12 x2 2 xl + 720x2x1 + 288 x2 2 - 12 xl 3 - 672 xl 2 

- 10560x1 + 1920x2 

"y[3,2]=", 145152x1 - 497664x2- 28 x2xl 3 + 1632 x2 2 xl - 528 x1 3 + 384 x2 3 
+ 72 x2 2 xl 2 - 8 xl 4 + 12 x2 3 xl - 2496x2x1 2 + x2 2 xl 3 - 73920x2x1 
+ 7680 x2 2 - 5760x1 2 
"y[3,3]=", 23553024x1 + 21829632x2+ 96 x2 3 xl 2 - 2160 x2xl 3 - 233280 x2 2 xl 
+ 62400x1 3 + 19200 x2 3 - 5760 x2 2 xl 2 + 480x1 4 - 24 x2x1 4 + 2880 x2 3 xl 

- 34560 x2 xl 2 - 48 x2 2 xl 3 + x2 3 xl 3 + 480 x2 4 + 1435392 x2 xl + 12 x2 4 xl 

- 2575872 x2 2 + 2345472x1 2 

> 
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